Simulation of N particle Interaction via Lennard Jones potential using python

Harsimran Singh

B.Sc-(Honours)-Physics-(III), Panjab University , Chandigarh

Abstract/Aim

Simulation of N particle interaction has been key for many advances in fields like Solid State Physics, Molecular Dynamics, Astrophysics and much more. As many times, systems under observation do not contain one atom or molecule but many of them arranged in random way where each configuartion has its own probability of occurence. In all such system one has to consider N body interaction inevitably. Due to such profound occurence, N body interaction simulation is strong tool which can be used in many fields like plasma physics. One use of N body interaction simulation had been in study of behaviour of gases and estimation of properties of atoms from behaviour of bulk. One of its example has been simulation of behaviour of argon which has been explained quite well by lennard jones potential. Due to its short range nature, It is the potential for which techniques of Molecular dynamics simulation works really well and it's computational cost is less. In this project, Simulation for Argon atoms interacting via l-J potential has been done using python3.

Introduction

Lennard Jones Potential

Lennard jones potential is short range potential which is attractive over large distancces but repulsive at very short distance. Form of this potential is

V = 4 ε ({(\frac{σ}{r})}^{12} - {(\frac{σ}{r})}^{6})

In this , minimum valuue of potential is at $r = 2^{1 / 6} σ$ and $V l j = - ε$

lennard Jones potential is the one of the most used potential while modeling for basic properties of neutral atoms and molecules, and both attractive and repulsive parts which account for van der waal forces and interatomic repulsion respectively . One of the most important advantage of using lj potential is because it is short range potential and dies out rapidly, which in turn opens the gateway for one simulation technique which makes computation times much lower than that of any other long range potential.

Lennard jones potential as function of r

N particle Interaction

In system containing more than one or two atoms or molecules, one has to consider interaction of each molecule with N-1 other atoms. From perspective of quantum mechanics, this is translated by Schroedinger wave equation.

i ℏ \frac{\partial ψ}{\partial t} = H ψ

here hamiltonion include potential due to all molecules and psi represent state of system. In classical mechanics, equation governing time evolution of system can be written as

\frac{d^{2} r_{i}}{d t^{2}} = \frac{1}{m_{i}} \sum_{i \neq j}^{j} f_{i j}

For simulation and computation we will be using classical equation and lj potential as governing equation of motion.

F = - \nabla V_{l j}

Method

Main aim of doing simulation is to compute position r and velocity v of each atom for later time t. this is done by computing position and velocity after each time dt and then repeating this process untill we reach desired state or time. Main part of program is to calculation of r(t+h) of each atom from previous configuation. Here we used two algorithm's. However only one of them has been used mostly used while other was used for confirmation of stability and uniqueness of solution.

Algorithm

Predictor Corrector method

Predictor corrector method is also known as modified euler method. It has been used excessively in domain of molecular dynamics because of its simplicity, reliability and lesser computation cost. This method had been used for simulating Argon behaviour (Rahman, et al, 1964 ). this algorithm is as:

Assuming we know position of an atom at time t=t-h , which we write as r ^n-1 , and position , velocity, accelaration at t=t, which is r ⁿ , v ⁿ, a ⁿ, now first we predict r ⁿ⁺¹ and then refine it later on as

r_{e}^{n + 1} = r^{n - 1} + 2 * △ t * v^{n}

from estimated position , we compute estimated accelaration a_e ⁿ⁺¹. After this, we calculate corrected value of r and v as

v^{n + 1} = v^{n} + \frac{1}{2} △ t * (a^{n} + a_{e}^{n + 1})

r^{n + 1} = r^{n} + \frac{1}{2} △ t * (v^{n} + v^{n + 1})

a ⁿ⁺¹can now be computed from calculated r ⁿ⁺¹. and process goes on.

Runge Kutte method

RK method is quite known method, But due to more computational cost it is only used here as tester for stability and uniquesness of solution provided by predictor corrector method.

k_{0} = △ t^{2} f (r_{n})

k 1 = △ t^{2} f (r_{n} + \frac{1}{2} △ t * v_{n} + \frac{1}{8} k_{0})

k 2 = △ t^{2} f (r_{n} + △ t * v_{n} + \frac{1}{2} k_{1})

r_{n + 1} = r_{n} + △ t * v_{n} + \frac{1}{6} (k_{0} + 2 k_{1})

v_{n + 1} = v_{n} + \frac{1}{6 △ t} (k_{0} + 4 k_{1} + k_{2})

Interaction¹

All pair method

It is the simplest to implement, but extremely inefficient when the interaction range r c is small compared with the linear size of the simulation region. All pairs of atoms must be examined because, owing to the continual rearrangement that characterizes the fluid state, it is not known in advance which atoms actually interact. Although testing whether atoms are separated by less than r c is only a part of the overall interaction computation, the fact that the amount of computation needed grows as O(N ²) rules out the method for all but the smallest values of N . Techniques for reducing this growth rate to a more acceptable O(N ) level, often used in tandem, will be discussed here. A schematic summary of the methods appears in Figure 2.

All pair method

Neighbour List method

In case of N particle interaction, computatiions are required for N(N-1)/2 pairs which results in computation varying roughly as N². For decreasing this dependence , In case of short range potential , assign neighbour list to all atoms , by which only those present in neighbours are considered for force calculation and this list can be updated after 10-20 cycles which makes computation time vary roughly as N and not as N².

figure depicting neighbor list method

Time saving techniques

To reduce time of computation which has been major drawbac of python3 , different compiler based python was used named pypy. It proved to make code run 6 times faster than on normal python, but was not capable of plotting and animation and thus all the data was once computed and saved and then was plottted and animate dusing python3.

Setting the parameters

for Argon, mass of an atom is 39.95x1.67x10^-27 Kg. In Lennard Jones Potential, we have $σ = 3.4 A n g s t o r m$ (Rahman, et al, 1964),(Verlet, et al, 1967).

In order to make compuatation easier following units have been adopted:

length-3.4 Amgstorm ( $σ$ )

Time- 10^-14 sec

$ε = K_{b} T_{a} (T_{a} = 120 K)$

from now on we will call our length unit X and time unit C (not to be confused by speed of light).

In these units, Boltzmann Constant is

K_boltzmann=1.064 x 10^-31 Kg X²/C²

Force equation in pair is

F = 24 ε \frac{1}{r^{2}} (\frac{2}{r^{12}} - \frac{1}{r^{6}}) ((x - x_{0}) i + (y - y_{0}) j + (z - z_{0}) k)

Thus This equation becomes:

\frac{d^{2} \vec{r}}{d t^{2}} = 24 \frac{ε}{M} \frac{1}{r^{2}} (\frac{2}{r^{12}} - \frac{1}{r^{6}}) ((x - x_{0}) i + (y - y_{0}) j + (z - z_{0}) k)

In our present units, for Argon, This eqquation becomes

\frac{d^{2} \vec{r}}{d t^{2}} = 0.0004593 \frac{1}{r^{2}} (\frac{2}{r^{12}} - \frac{1}{r^{6}}) ((x - x_{0}) i + (y - y_{0}) j + (z - z_{0}) k)

Law of Conservation Of Energy as test of Reliabilty

In Simulation, One is not aware of final output and depends entirely on simulation itself for validation. In such cases, Whether our Algorithm has errors under control or not is hard to determine in exact way because one do not know theoretical values of such system. In order to accomplish this, Method of Law of Conservation of Energy is used. System Initial total energy is compared with system final Energy to get the Idea about Uniqueness and stability of Code.

Coding Essentials

Code has been made using Object Oriented Programming Paradigm. It is due to the fact that we need to take care of properties of particles and which are huge in number , thus makig it quite redundant to implement in procedural programming. As, Python is object oriented language along with its ability to loop over objects, It became primary choice. Some basics are discussed here.

Object Oriented Programming

Object oriented programming allows us to make objects with specific attributes and methods, which we used to create particles. As our particles have properties like position, velocity, accelaration associated making a class particle was inevitabe part. Here we discuss brief method of making class in python3.

                                       Class Particle():   ///// Starting Class Definition
    def __init__(self,position,velocity):

Here position and velocity are additional Arguments, which one hould enter while creating an object.

                                       class particle():
    def __init__(self, position,velocity):
        self.position=position   /// Defining attributes of object
        self.velocity=velocity
        self.accelaration=(0,0,0)
         
                               

Here self.X is attribute with name X which belongs to the object , in our case particle. As it is evident that every Object attribute need not to be in definition of __init__().

Now we will se how to add function in our class and objects. Here Two distinction occurs. In code two kind of Methods(Functions in Class are known as methods).

1. Methods which belong to objects. That is every object have access to them. They are defined as normal functions.

                                       def method_name(self,*Args): // self is our object. , *Args are other user defined arguments
///Using This method
object.method_name(*Args)
         
                               

2. Methods which belong to class only and cannot be accessed by objects are Classmethods. They are also defined in same way except we have to specify that it is class method

                                       @classmethod
def classmethod_name(cls, *Args):   \\\ cls is argument used to denote class, *Args are other user defined arguments.
\\\Using class Method. In our case class name is Particle.
Particle.classmethod_name(*Args)
         
                               

For More Info about Python classes and objects. See-(Learn Python-The Hard way)

Code

we first define class particle which was used to create particle objects with attributes required in computation for each particle ,object and class methods.

Defining Class

                                       ### Defining class particle for creating particle objects
class particle():
    ### List to store handle of each particle created
    objects = []
    
    # time variable for keeping track of time of simulation
    time=0
    
    # defining initiating function of class, which creates object
    def __init__(self, position, velocity):
        
        #position, velocity , accelaration and neighbour list attribute of particle
        self.position = position
        self.velocity = velocity
        self.accelaration = (0, 0, 0)
        self.neighbours = []
        
        
        # temporary position and velocity variable used in algorithm
        self.position1=position
        self.velocity1=velocity
        self.velocity2=velocity
        self.velocity3=velocity
        self.accelaration1=(0,0,0) 
        self.da=(0,0,0)
        
        # adding object handle in objects list
        particle.objects.append(self)
    # method for updating position and velocity by 1 step
    
    @classmethod
    def update(cls):
        # step size
        t=0.1
        
        # half of step size required for estimation
        th=t/2
        
        # updating neighbours
        if particle.time%100==0:
            particle.update_neighbours()
            
        # updating accelaration before updating coordinates
        particle.update_ac()
        
        # calculating position and velocity
        
        for obj in cls.objects:
            # position at n-1 step
            xn1, yn1, zn1 = obj.position1
            
            # velocity
            vx,vy,vz=obj.velocity
            
            # creating box by condition
            if xn1>=2.7:
                vx,vy,vz=(-abs(vx),vy,vz)
              
            if yn1>=2.7:
                vx,vy,vz=(vx,-abs(vy),vz)         
               
            if zn1>=2.7:
                vx,vy,vz=(vx,vy,-abs(vz))
              
            if xn1<=-2.7:
                vx,vy,vz=(abs(vx),vy,vz)
       
            if yn1<=-2.7:
                vx,vy,vz=(vx,abs(vy),vz)
          
            if zn1<=-2.7:
                vx,vy,vz=(vx,vy,abs(vz))
                
            
            obj.velocity=(vx,vy,vz)
            
            # storing current velocity and accelaration and position in temporary variable
            obj.velocity1=obj.velocity
            obj.accelaration1=obj.accelaration
            obj.position1=obj.position
            
            # updating poition for estimated value at n+1 step
            obj.position = (xn1 + 2*t * vx, yn1 + 2*t * vy, zn1 + 2*t * vz)
           
        # updating accelaration for estimated location
        particle.update_ac()
        
        # correcting estimate 
        for obj in particle.objects:
            vxn,vyn,vzn=obj.velocity
            
            obj.velocity=(vxn+(obj.accelaration[0]+obj.accelaration1[0])*th,vyn+(obj.accelaration[1]+obj.accelaration1[1])*th,vzn+(obj.accelaration[2]+obj.accelaration1[2])*th)  
            obj.position=(obj.position1[0]+(vxn+obj.velocity[0])*th,obj.position1[1]+(vyn+obj.velocity[1])*th,obj.position1[2]+(vzn+obj.velocity[2])*th)
            
        # incrementing time 
        particle.time+=1
        
    # defining force function
    def force(self):
        ep = 0.0004593 # Well depth kb_m*T , T= 120K
        d = 1 # sigma
        f = (0, 0, 0)
        x0, y0, z0 = self.position
        for obj2 in self.neighbours:
            x, y, z = obj2.position
            r=((x-x0)**2+(y-y0)**2 +(z-z0)**2)
            if r<7.84:
            
                p = ep * (2 * ((1)/(r**7)) - 1 * ((1)/(r**4)))
                f = (f[0] -p * (x - x0), f[1] - p * (y - y0), f[2] - p * (z - z0))
       
        # changing accelaration so as to not calculate force for future config.
        self.accelaration=f
        
        
    # creating update neighbour function
    @classmethod
    def update_neighbours(cls):
        R_critical=3 # 2.5*Sigma
        
        for obj in cls.objects:
            obj.neighbours=[]
            x,y,z=obj.position
            for obj2 in cls.objects:
                x1,y1,z1=obj2.position
                r=((x-x1)**2+(y-y1)**2+(z-z1)**2)**0.5
                if r<R_critical and r>0:
                    obj.neighbours.append(obj2)
    @classmethod
    def update_neighbours1(cls):
        R_critical=2.8 # 2.5*Sigma
        
        for obj in cls.objects:
            obj.neighbours=[]
            x,y,z=obj.position
            for obj2 in cls.objects:
                x1,y1,z1=obj2.position
                r=((x-x1)**2+(y-y1)**2+(z-z1)**2)**0.5
                if r<R_critical and r>0:
                    obj.neighbours.append(obj2)
    # calculating energy
    @classmethod
    def energy(cls):
        e = 0
        for obj in cls.objects:
            vx, vy, vz = obj.velocity
            e = e + vx**2 + vy**2 + vz**2
        return e
        
    # creating function to update acceleration
    @classmethod
    def update_ac(cls):
        for obj in cls.objects:
            obj.force()
    
    x_list = []
    y_list = []
    z_list = []
    
    # function to get position of all particles
    @classmethod
    def getcoordinate(cls):
        cls.x_list = []
        cls.y_list = []
        cls.z_list = []
        for obj in cls.objects:
            x, y, z = obj.position
            cls.x_list.append(x)
            cls.y_list.append(y)
            cls.z_list.append(z)
        return cls.x_list, cls.y_list, cls.z_list
        
    # function to get velocities of all particles
    @classmethod
    def getvelocity(cls):
        cls.x_list = []
        cls.y_list = []
        cls.z_list = []
        for obj in cls.objects:
            x, y, z = obj.velocity
            cls.x_list.append(x)
            cls.y_list.append(y)
            cls.z_list.append(z)
        return cls.x_list, cls.y_list, cls.z_list
        
    # functions for generating velocity distribution data initial and final   
    @classmethod
    def velocity_distribution(cls):
        v=[]
        for obj in cls.objects:
            v.append((obj.velocity[0]**2+obj.velocity[1]**2+obj.velocity[2]**2)**0.5)
        np.save('velocity_data.npy',v,allow_pickle=True)
    @classmethod
    def velocity_distribution1(cls):
        v1=[]
        for obj in cls.objects:
            v1.append((obj.velocity[0]**2+obj.velocity[1]**2+obj.velocity[2]**2)**0.5)
        np.save('velocity_data1.npy',v1,allow_pickle=True)
    
    
    # function for calculating potential energy due to neighbours
    @classmethod
    def potential(cls):
        potent=0
        for obj in cls.objects:
            x0,y0,z0= obj.position
            for obj2 in obj.neighbours:
                x, y, z = obj2.position
                r=((x-x0)**2+(y-y0)**2 +(z-z0)**2)**(0.5)
                potent+=480*(1/(r**12)  -  1/(r**6))
        potent=potent/2
        return potent
     
    # function to calculate potential energy due to all particles
    @classmethod
    def potential1(cls):
        potent=0
        for obj in cls.objects:
            x0,y0,z0= obj.position
            for obj2 in cls.objects:
                x, y, z = obj2.position
                r=((x-x0)**2+(y-y0)**2 +(z-z0)**2)**(0.5)
                if r>0:
                    potent+=480*(1/(r**12)  -  1/(r**6))
        return potent/2
         
                      defining class         

In order to create particles with given position and random velocity, create function was created as follows:

                                       def create():
    lists = []
    v_da1=[]
    x_da1=[]
    for i in range(125):
        
        x = 0.001*np.random.randint(-1000,1000)
        y = 0.001*np.random.randint(-1000,1000)
        z = 0.001*np.random.randint(-1000,1000)
        
        vx = 0.0001*rd.randint(-100,100)*0.867
        vy = 0.0001*rd.randint(-100,100)*0.867
        vz = 0.0001*rd.randint(-100,100)*0.867
        x_da1.append((x,y,z))
        v_da1.append((vx,vy,vz))
        lists.append(particle((x, y, z), (vx, vy, vz)))
    np.save("v_da1.npy",v_da1)
    np.save("x_da1.npy",x_da1)
         
                      Create function         

In order to calculate local density, whole volume was divided into 64 parts. and number of particle inside that box were counted and used in calculating standard deviation of number of particles in each box. and this deviation was calculated for every time.

                                       import math
density=np.zeros((4,4,4))
density_data=[]
for data in position:
    density=np.zeros((4,4,4))
    x,y,z=data
    for number in range(len(x)):
        density[1+math.ceil(x[number]/7.51),1+math.ceil(y[number]/7.51),1+math.ceil(z[number]/7.51)]+=1
    density_data.append(np.std(density)/np.mean(density))
np.save('density.npy',density_data,allow_pickle=True)
         
                      Denisty Deviation         

At the end we have animation and plotting code

                                       import os
import numpy as np
print('\n Plotting and animation Begins\n\n')
os.chdir('')
a=np.load('data_set.npy',allow_pickle=True)
v1=np.load('velocity_data1.npy',allow_pickle=True)
v=np.load('velocity_data.npy',allow_pickle=True)
T_data=np.load('vel_data.npy',allow_pickle=True)
step=0.0005
final=0.02
mean_v=np.mean(v)
mean_v1=np.mean(v1)
meansv=np.mean(v**2)**0.5
meansv1=np.mean(v1**2)**0.5
M=39.95*1.67*10**(-27)
Kb=1.064*10**(-32)
#Temperature
Ti=(M/(3*Kb))*(np.mean(v1**2))
Tf=(M/(3*Kb))*(np.mean(v**2))
vr=np.arange(0,final,step)
Maxw=125*(4*np.pi*vr**2*step)*((M/(2*np.pi*Kb*Tf))**(1.5))*np.exp(-(M/(2*Kb*Tf))*vr**2)
print('Mean v:',np.mean(v))
import matplotlib.pyplot as plt
plt.plot(T_data)
plt.title('Variation Of temperature with Time')
plt.ylabel('Temparature')
plt.ylim(60,1.1*np.max(T_data))
plt.figure()
plt.subplot(121)
a1,b1,c1=plt.hist(v1,bins=np.arange(0,final,step))
ind=np.argmax(a1)
plt.title("Initial Velocity Distribution")
plt.subplot(122)
a2,b2,c2=plt.hist(v,bins=np.arange(0,final+step,step))
ind2=np.argmax(a2)
if np.max(a1)>np.max(a2):
    limit=np.max(a1)
else:
    limit=np.max(a2)
plt.title("Final Velocity Distribution")
plt.text(0.012,1.4*limit-4,"Temperature: {:.2f}".format(Tf), bbox=dict(boxstyle='round',facecolor='white'))
plt.text(0.012,1.4*limit-5,"Mean v : {:.4f}".format(mean_v)+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
plt.text(0.012,1.4*limit-6,"RMS v : {:.4f}".format(meansv)+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
plt.text(0.012,1.4*limit-7,"MP v : {:.4f}-{:.4f}".format(b2[ind2],b2[ind2+1])+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
plt.plot(vr,Maxw)
plt.legend(["Maxwell Distribution","Observed Distribution"])
plt.subplot(121)
plt.ylim(0,1.4*limit)
plt.text(0.012,1.4*limit-4,"Temperature: {:.2f}".format(Ti), bbox=dict(boxstyle='round',facecolor='white'))
plt.text(0.012,1.4*limit-5,"Mean v : {:.4f}".format(mean_v1)+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
plt.text(0.012,1.4*limit-6,"RMS v : {:.4f}".format(meansv1)+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
plt.text(0.012,1.4*limit-7,"MP v : {:.4f}-{:.4f}".format(b1[ind],b1[ind+1])+"X/c", bbox=dict(boxstyle='round',facecolor='white'))
plt.subplot(122)
plt.ylim(0,1.4*limit)
plt.show()
print('###################')
density=np.load('density.npy',allow_pickle=True)
plt.plot(density)
plt.title('Variation of Deviation in Local Density with time') 
plt.ylim(0.1,1.1*np.max(density))
plt.show()
print('wait a second')
def update_coordinate(i, scatter):
    x, y, z = a[i]
    ax.clear()
    ax.set_xlim3d([-3,3])
    ax.set_ylim3d([-3,3])
    ax.set_zlim3d([-3,3])
    ax.set_title('Interaction Animation')
    scatter[0]=ax.scatter3D(x,y,z,color='green',marker=".")
import mpl_toolkits.mplot3d.axes3d as mp
import matplotlib.animation as animation
Writer = animation.writers['ffmpeg']
writer = Writer(fps=10, metadata=dict(artist='Me'), bitrate=2800)
fig=plt.figure()
ax = plt.axes(projection ="3d")
ax.set_xlim3d([-3,3])
ax.set_ylim3d([-3,3])
ax.set_zlim3d([-3,3])
ax.set_title("Scatter Animation")
x,y,z=a[0]
scatter=[ax.scatter3D(x,y,z,color='green')]
anim=animation.FuncAnimation(fig,update_coordinate,len(a),fargs=(scatter,),interval=20)
anim.save("LJ_Potential.mp4",writer=writer)
plt.show()
         
                      Animation and plot code         

one script file was written to run all codes on different platform as required.

                                       pypy3 main.py
pypy3 density_distribution.py
python3 animate.py
         
                      Script file for terminal         

Result

Simulation was carried out for two different Densities of Argon p=1.3954 g/cc and p=0.6977 g/cc. foe three different Temperatures each.

Initial T=79 K , $ρ = 1.3954 g / c c$

Initial Temperature=79 K

Step Size=0.1

Integration Time: 50,000 x 0.1

Time taken for execution : 205.84 sec

Initial average Kinetic Energy=118.67 Kb(Boltzmann Constant)

Initial average potential energy of an atom = -221 Kb

Final Kinetic Energy= 133.88 Kb

Final average potential energy of an atom = -235.44 Kb

Total error in Temperature = 0.5153 K

Total Error in energy in Kb = 0.773

Animation

Velocity Distribution

Variation Of Temperature with time

Variation of local density deviation with time

Results of simulation

Initial T=163 K , $ρ = 1.3954 g / c c$

Initial Temperature=163.3 K

Step Size=0.1

Integration Time: 50,000 x 0.1

Time taken for execution : 50.95 sec

Initial average Kinetic Energy=244.95 Kb(Boltzmann Constant)

Initial average potential energy of an atom = -221 Kb

Final Kinetic Energy= 136.82 Kb

Final average potential energy of an atom = -106.76 Kb

Total error in Temperature = 4.073 K

Total Error in energy in Kb = 6.1105

Animation

Velocity Distribution

Variation of Temperature with Time

Variation of Local density deviation with time

Initial T=202 K, $ρ = 1.3954 g / c c$

Initial Temperature=202.3 K

Step Size=0.1

Integration Time: 50,000 x 0.1

Time taken for execution : 29.19 sec

Initial average Kinetic Energy=303.5 Kb(Boltzmann Constant)

Initial average potential energy of an atom = -221 Kb

Final Kinetic Energy= 140.91 Kb

Final average potential energy of an atom = -49.62 Kb

Total error in Temperature = 5.86 K

Total Error in energy in Kb = 8.739

Animation

Velocity Distribution

Variation of Temperature with Time

Variation of Local density deviation with time

Initial T=98 K, $ρ = 0.697 g / c c$

Initial Temperature=97.6 K

Step Size=0.1

Integration Time: 50,000 x 0.1

Time taken for execution : 19.16 sec

Initial average Kinetic Energy=146.41 Kb(Boltzmann Constant)

Initial average potential energy of an atom = -60.37 Kb

Final Kinetic Energy= 117.93 Kb

Final average potential energy of an atom = -26.16 Kb

Total error in Temperature = 3.81 K

Total Error in energy in Kb = 5.72

Animation

Velocity Distribution

Variation of Temperature with Time

Variation of Local density deviation with time

Initial T=168 , $ρ = 0.697 g / c c$

Initial Temperature=168.3 K

Step Size=0.1

Integration Time: 50,000 x 0.1

Time taken for execution : 14.60 sec

Initial average Kinetic Energy=252.8 Kb(Boltzmann Constant)

Initial average potential energy of an atom = -60.37 Kb

Final Kinetic Energy= 211.75 Kb

Final average potential energy of an atom = -13.41 Kb

Total error in Temperature = 3.93 K

Total Error in energy in Kb = 5.91

Velocity Distribution

Variation of Temperature with Time

Variation of Local density deviation with time

Initial T=214 K , $ρ = 0.697 g / c c$

Initial Temperature=214.29 K

Step Size=0.1

Integration Time: 50,000 x 0.1

Time taken for execution : 13.34 sec

Initial average Kinetic Energy=321.44 Kb(Boltzmann Constant)

Initial average potential energy of an atom = -60.37 Kb

Final Kinetic Energy= 280.70 Kb

Final average potential energy of an atom = -13.84 Kb

Total error in Temperature = 3.87 K

Total Error in energy in Kb = 5.80

Animation

Velocity Distribution

Variation of Temperature with Time

Variation of Local density deviation with time

References

Reference

1

Context
Rahman, A. (1964) “Correlations in the Motion of Atoms in Liquid Argon.”

1
Verlet, Loup (1967) “Computer ‘Experiments’ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules.”

1

Source

Fig 1a: https://medium.com/search?q=design
Fig 2: D.C. Rapport, Cambridge University Press, The Art of molecular Dynamics Simulation
Fig 3: The Art of Molecular Dynamics by D.C. Rapport

Manage ContentSee Navigation

0
0 Simulation of N particle Interaction via Lennard Jones potential using python
0 Abstract/Aim
0 1 Introduction
0 1.1 Lennard Jones Potential
0 1.2 N particle Interaction
0 2 Method
0 2.1 Algorithm
0 2.1.1 Predictor Corrector method
0 2.1.2 Runge Kutte method
0 2.2 Interaction1
0 2.2.1 All pair method
0 2.2.2 Neighbour List method
0 2.3 Time saving techniques
0 2.4 Setting the parameters
0 2.5 Law of Conservation Of Energy as test of Reliabilty
0 2.6 Coding Essentials
0 2.6.1 Object Oriented Programming
0 2.7 Code
0 3 Result
0 3.1 Initial T=79 K , $\rho=1.3954\;g/cc$ <math xmlns="http://www.w3.org/1998/Math/MathML"><mi>ρ</mi><mo>=</mo><mn>1</mn><mo>.</mo><mn>3954</mn><mo> </mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=1.3954\;g/cc
0 3.2 Initial T=163 K , $\rho=1.3954\;g/cc$ <math xmlns="http://www.w3.org/1998/Math/MathML"><mi>ρ</mi><mo>=</mo><mn>1</mn><mo>.</mo><mn>3954</mn><mo> </mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=1.3954\;g/cc
0 3.3 Initial T=202 K, $\rho=1.3954\;g/cc$ <math xmlns="http://www.w3.org/1998/Math/MathML"><mi>ρ</mi><mo>=</mo><mn>1</mn><mo>.</mo><mn>3954</mn><mo> </mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=1.3954\;g/cc
0 3.4 Initial T=98 K, $\rho=0.697\;g/cc$ <math xmlns="http://www.w3.org/1998/Math/MathML"><mi>ρ</mi><mo>=</mo><mn>0</mn><mo>.</mo><mn>697</mn><mo> </mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=0.697\;g/cc
0 3.5 Initial T=168 , $\rho=0.697\;g/cc$ <math xmlns="http://www.w3.org/1998/Math/MathML"><mi>ρ</mi><mo>=</mo><mn>0</mn><mo>.</mo><mn>697</mn><mo> </mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=0.697\;g/cc
0 3.6 Initial T=214 K , $\rho=0.697\;g/cc$ <math xmlns="http://www.w3.org/1998/Math/MathML"><mi>ρ</mi><mo>=</mo><mn>0</mn><mo>.</mo><mn>697</mn><mo> </mo><mi>g</mi><mo>/</mo><mi>c</mi><mi>c</mi></math>\rho=0.697\;g/cc

0
0 Fig 1 : Lennard jones potential as function of r
0 Fig 2 : All pair method
0 Fig 3 : figure depicting neighbor list method
0 Video 1 : Animation
0 Fig 4 a : Velocity Distribution
0 Fig 4 b : Variation Of Temperature with time
0 Fig 4 c : Variation of local density deviation with time
0 Video 2 : Animation
0 Fig 5 a : Velocity Distribution
0 Fig 5 b : Variation of Temperature with Time
0 Fig 5 c : Variation of Local density deviation with time
0 Video 3 : Animation
0 Fig 6 a : Velocity Distribution
0 Fig 6 b : Variation of Temperature with Time
0 Fig 6 c : Variation of Local density deviation with time
0 Video 4 : Animation
0 Fig 7 a : Velocity Distribution
0 Fig 7 b : Variation of Temperature with Time
0 Fig 7 c : Variation of Local density deviation with time
0 Video 5 :
0 Fig 8 a : Velocity Distribution
0 Fig 8 b : Variation of Temperature with Time
0 Fig 8 c : Variation of Local density deviation with time
0 Video 6 : Animation
0 Fig 9 a : Velocity Distribution
0 Fig 9 b : Variation of Temperature with Time
0 Fig 9 c : Variation of Local density deviation with time

References

Hyperlink

Notes

Glossary

Source

J Hill, C Bargmann

Cited in this article

Section Title

Cited in this article

Equation #

Cited in this article

Figure Caption

Cited in this article

Table Caption

Cited in this article

Section Title

Cited in this article

Section Title

Cited in this article

Chat